Technical details
library(GeoPressureR)
library(leaflet)
library(leaflet.extras)
library(raster)
library(dplyr)
library(ggplot2)
library(kableExtra)
library(plotly)
library(GeoLocTools)
setupGeolocation()
knitr::opts_chunk$set(echo = FALSE)
load(paste0("../data/1_pressure/", params$gdl_id, "_pressure_prob.Rdata"))
load(paste0("../data/2_light/", params$gdl_id, "_light_prob.Rdata"))
load(paste0("../data/3_static/", params$gdl_id, "_static_prob.Rdata"))
load(paste0("../data/4_basic_graph/", params$gdl_id, "_basic_graph.Rdata"))
All the results produced here are generated with (1) the raw geolocator data, (2) the labeled files of pressure and light and (3) the parameters listed below.
kable(gpr) %>% scroll_box(width = "100%")
| gdl_id | crop_start | crop_end | thr_dur | extent_N | extent_W | extent_S | extent_E | map_scale | map_max_sample | map_margin | prob_map_s | prob_map_s_calib | prob_map_thr | shift_k | kernel_adjust | calib_lon | calib_lat | calib_1_start | calib_1_end | calib_2_start | calib_2_end | calib_2_lon | calib_2_lat | prob_light_w | thr_prob_percentile | thr_gs | RingNo | scientific_name | common_name | mass | wing_span | Color |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 16JL | 2016-07-01 | 2016-11-12 20:00:00 | 24 | 51 | -18 | 4 | 16 | 5 | 300 | 30 | 1 | 2 | 0.9 | 0 | 1.4 | 8.68315 | 46.54263 | 2016-07-01 | 2016-09-18 | NA | NA | NA | NA | 0.1 | 0.9 | 120 | N566642 | Oenanthe oenanthe | Northern wheatear | NA | NA | NA |
The labeling of pressure data is illustrated with this figure. The black dots indicates the pressure datapoint not considered in the matching. Each stationary period is illustrated by a different colored line.
pressure_na <- pam$pressure %>%
mutate(obs = ifelse(isoutliar | sta_id == 0, NA, obs))
p <- ggplot() +
geom_line(data = pam$pressure, aes(x = date, y = obs), colour = "grey") +
geom_point(data = subset(pam$pressure, isoutliar), aes(x = date, y = obs), colour = "black") +
# geom_line(data = pressure_na, aes(x = date, y = obs, color = factor(sta_id)), size = 0.5) +
geom_line(data = do.call("rbind", shortest_path_timeserie) %>% filter(sta_id > 0), aes(x = date, y = pressure0, col = factor(sta_id))) +
theme_bw() +
scale_colour_manual(values = col) +
scale_y_continuous(name = "Pressure(hPa)")
ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)
raw_geolight <- pam$light %>%
transmute(
Date = date,
Light = obs
)
lightImage(tagdata = raw_geolight, offset = 0)
tsimagePoints(twl$twilight,
offset = 0, pch = 16, cex = 1.2,
col = ifelse(twl$deleted, "grey20", ifelse(twl$rise, "firebrick", "cornflowerblue"))
)
abline(v = gpr$calib_2_start, lty = 1, col = "firebrick", lwd = 1.5)
abline(v = gpr$calib_1_start, lty = 1, col = "firebrick", lwd = 1.5)
abline(v = gpr$calib_2_end, lty = 2, col = "firebrick", lwd = 1.5)
abline(v = gpr$calib_1_end, lty = 2, col = "firebrick", lwd = 1.5)
The probability map resulting from light data alone can be seen below.
li_s <- list()
l <- leaflet(width = "100%") %>%
addProviderTiles(providers$Stamen.TerrainBackground) %>%
addFullscreenControl()
for (i_r in seq_len(length(light_prob))) {
i_s <- metadata(light_prob[[i_r]])$sta_id
info <- pam$sta[pam$sta$sta_id == i_s, ]
info_str <- paste0(i_s, " | ", info$start, "->", info$end)
li_s <- append(li_s, info_str)
l <- l %>% addRasterImage(light_prob[[i_r]], opacity = 0.8, colors = "OrRd", group = info_str)
}
l %>%
addCircles(lng = gpr$calib_lon, lat = gpr$calib_lat, color = "black", opacity = 1) %>%
addLayersControl(
overlayGroups = li_s,
options = layersControlOptions(collapsed = FALSE)
) %>%
hideGroup(tail(li_s, length(li_s) - 1))
We can compare light and pressure location at long stationary stopover (>5 days). By assuming the best match of the pressure to be the truth, we can plot the histogram of the zenith angle and compare to the fit of kernel density at the calibration site.
raw_geolight <- pam$light %>%
transmute(
Date = date,
Light = obs
)
dur <- unlist(lapply(static_prob_marginal, function(x) difftime(metadata(x)$temporal_extent[2],metadata(x)$temporal_extent[1], units = "days" )))
long_id <- which(dur>5 & seq(1,length(dur))<=length(shortest_path_timeserie))
par(mfrow = c(2, 3))
for (i_s in long_id){
twl_fl <- twl %>%
filter(!deleted) %>%
filter(twilight>shortest_path_timeserie[[i_s]]$date[1] & twilight<tail(shortest_path_timeserie[[i_s]]$date,1))
sun <- solar(twl_fl$twilight)
z_i <- refracted(zenith(sun, shortest_path_timeserie[[i_s]]$lon[1], shortest_path_timeserie[[i_s]]$lat[1]))
hist(z_i, freq = F, main = paste0("sta_id=",i_s, " | ",nrow(twl_fl),"twls"))
lines(fit_z, col = "red")
xlab("Zenith angle")
}
Similarly, we can plot the line of sunrise/sunset at the best match of pressure (yellow line) and compare to the raw and labeled light data.
lightImage(
tagdata = raw_geolight,
offset = gpr$shift_k / 60 / 60
)
tsimagePoints(twl$twilight,
offset = gpr$shift_k / 60 / 60, pch = 16, cex = 1.2,
col = ifelse(twl$deleted, "grey20", ifelse(twl$rise, "firebrick", "cornflowerblue"))
)
for (ts in shortest_path_timeserie){
if (!is.null(ts)){
twl_fl <- twl %>%
filter(twilight>ts$date[1] & twilight<tail(ts$date,1))
if (nrow(twl_fl)>0){
tsimageDeploymentLines(twl_fl$twilight,
lon = ts$lon[1], ts$lat[1],
offset = gpr$shift_k / 60 / 60, lwd = 3,col = adjustcolor("orange", alpha.f = 0.5))
}
}
}
To visualize the path on GeoPressureViz, you will need to also load the pressure and light probability map and align them first with the code below.
sta_marginal <- unlist(lapply(static_prob_marginal, function(x) raster::metadata(x)$sta_id))
sta_pres <- unlist(lapply(pressure_prob, function(x) raster::metadata(x)$sta_id))
sta_light <- unlist(lapply(light_prob, function(x) raster::metadata(x)$sta_id))
pressure_prob <- pressure_prob[sta_pres %in% sta_marginal]
light_prob <- light_prob[sta_light %in% sta_marginal]
The code below will open with the shortest path computed with the graph approach.
geopressureviz <- list(
pam_data = pam,
static_prob = static_prob,
static_prob_marginal = static_prob_marginal,
pressure_prob = pressure_prob,
light_prob = light_prob,
pressure_timeserie = shortest_path_timeserie
)
save(geopressureviz, file = "~/geopressureviz.RData")
shiny::runApp(system.file("geopressureviz", package = "GeoPressureR"),
launch.browser = getOption("browser")
)
| start | end | sta_id |
|---|---|---|
| 2016-07-01 00:00:00 | 2016-09-17 18:00:00 | 1 |
| 2016-09-17 20:00:00 | 2016-09-18 18:00:00 | 2 |
| 2016-09-19 03:30:00 | 2016-09-19 18:00:00 | 3 |
| 2016-09-20 02:30:00 | 2016-09-20 19:30:00 | 4 |
| 2016-09-20 21:00:00 | 2016-10-02 17:30:00 | 5 |
| 2016-10-03 00:30:00 | 2016-10-03 17:30:00 | 6 |
| 2016-10-03 23:00:00 | 2016-10-04 18:00:00 | 7 |
| 2016-10-05 04:30:00 | 2016-10-05 20:00:00 | 8 |
| 2016-10-06 05:00:00 | 2016-10-06 19:30:00 | 9 |
| 2016-10-06 20:30:00 | 2016-10-15 17:30:00 | 10 |
| 2016-10-16 05:30:00 | 2016-10-16 17:00:00 | 11 |
| 2016-10-16 22:00:00 | 2016-10-16 23:30:00 | 12 |
| 2016-10-17 06:00:00 | 2016-10-17 17:00:00 | 13 |
| 2016-10-18 05:00:00 | 2016-10-18 18:30:00 | 14 |
| 2016-10-18 22:30:00 | 2016-10-28 18:00:00 | 15 |
| 2016-10-29 05:30:00 | 2016-10-29 18:30:00 | 16 |
| 2016-10-30 02:00:00 | 2016-11-12 12:00:00 | 17 |